should you optimize?

software engineers should worry about other issues (such as good algorithm design and good implementations of those algorithms) before they worry about micro-optimizations” (Hyde 2009)

Do you have other sh*t to do? No? Please contact me and I’ll help you with that. Yes? You are among the 90% of R users whose first priority is not computer programming. The time spent optimizing code is often longer than the computing time actually saved. Use the simple solutions and you can get on with your research/life. (Ross 2013)

• but faster code enables interactivity, exploration, resampling-based methods …

categories of optimization

• low-hanging fruit
• parallelization
• (memoization)
• automatic differentiation
• faster languages

measuring code performance

• measure speed before you try to fix it
• complexity theory/big-O notation (e.g., \(f(x) = {\cal O}(x^{3/2})\) is useful if
  • you have some idea of the scaling constant
  • you plan to scale to large problems (or someone does)
• scaling dimensions: number of observations, parameters, clusters, dimensions …
• Amdahl’s Law: speeding up bottlenecks is what matters

benchmarking

• examples should be large enough not to be swamped by overhead
• … but small enough to run replicates and average
• system.time() for quick estimates
• rbenchmark, microbenchmark packages
• estimating scaling: theoretical or empirical (log-log plot)

scaling example

Modified from Brooks et al. (2017):

profiling

Example random-walk code (Ihaka, Temple Lang, and McArdle 2009):

rw2d2 <- function(n) {
    steps <- sample(c(-1, 1), n - 1, replace = TRUE)
    xdir <- sample(c(TRUE, FALSE), n - 1, replace = TRUE)
    xpos <- c(0, cumsum(ifelse(xdir, steps, 0)))
    ypos <- c(0, cumsum(ifelse(xdir, 0, steps)))
    return(list(x = xpos, y = ypos))
}

profiling run

• run multiple times to collect enough data/average over variation

Rprof("Rprof.out") ## start profiling
for (i in 1:100) {
    pos <- rw2d2(1e5)
}
Rprof(NULL) ## stop profiling

profiling results

source("https://raw.githubusercontent.com/noamross/noamtools/master/R/proftable.R")
proftable("Rprof.out", lines = 5)

 PctTime Call               
 38.18   sample > sample.int
 29.09   ifelse > which     
 18.18   ifelse             
  7.27   sample             
  3.64   c                  

Parent Call: rw2d2 > ...

Total Time: 1.1 seconds
Percent of run time represented: 96.4 %

profiling

• see also summaryRprof() (base-R), profvis package (RStudio fanciness)
• easiest to interpret if code is organized into functions
• profiling C++ code called from R is harder but not impossible

don’t write bad code

calculating \(\pi\) by Monte Carlo

bad code

non-vectorized, growing objects

piR_slow <- function(N) {
    res <- numeric(0)
    for (i in 1:N) {
        res <- c(res, as.numeric(runif(1)^2 + runif(1)^2 < 1.0))
    }
    4*sum(res)/N
}

better code

piR <- function(N) {
    x <- runif(N)
    y <- runif(N)
    d <- sqrt(x^2 + y^2)
    return(4 * sum(d < 1.0) / N)
}

benchmarking

rbenchmark::benchmark(
                bad = piR_slow(1e4),
                better = piR(1e4),
                columns = c("test", "replications", "elapsed", "relative"))

    test replications elapsed relative
1    bad          100  10.951  353.258
2 better          100   0.031    1.000

use good packages

• especially for data handling, I/O
  • data.table > tidyverse > base R
  • data formats: arrow, vroom, direct database access (dbplyr)
• collapse, xts, Rfast

sum-by-groups benchmark

• From https://h2oai.github.io/db-benchmark/
• 100M rows x 9 cols = 5 GB

get a bigger computer

parallelization types

• distributed: “embarrassingly parallel”
• multicore: multiple R processes
  • each one copies all objects
• threading: lightweight, shared-memory
• see also: parallel and HPC task view

distributed computing

• start individual runs, or chunks (e.g. simulations, bootstraps) as completely separate R jobs
• Compute Canada (or AWS or a big local workstation)
• jobs assigned by batch scheduler: High Performance Computing
• some tools within R: slurm, batchtools, futures.batchtools package

multicore

• spawn multiple jobs on a single machine (multiple cores)
• parallel package, foreach/doParallel

threading

• not available directly through R
• OpenMP
• parallel BLAS
  • base R comes with a robust but slow linear algebra library
  • can replace with optimized/parallel versions (command-line bullshittery required (Browne 2021))

BLAS benchmarks

• from https://csantill.github.io/RPerformanceWBLAS/:

AD basics

• automated chain rule (Griewank 1989; Griewank and Walther 2003)
• much faster and more robust than finite differences; faster than symbolic derivatives
• computing gradient takes \(\le 5 \times\) the effort of computing the objective function, regardless of dimension
• reverse-mode autodiff for \(f: {\mathbf R}^n \to {\mathbf R}\)
• widely used in modern ML (neural net back-propagation)
• toolboxes: TensorFlow (PyTorch, RTorch (Keydana 2023)), Stan (Chau 2022), (R)TMB, Julia tools …

a side note on optim()

• if you use a gradient-based optimizer (BFGS, L-BFGS-B) in optim() R will implement finite-difference gradients by default
• this is terrible (slow and fragile), although fine for simple problems
• … might as well use gradient-free methods (e.g. Nelder-Mead) for robustness
• providing gradients (somehow) can make a huge difference for tough optimization problems

AD example

deriv(~cos(log(1+x^2)), "x")

expression({
    .expr2 <- 1 + x^2
    .expr3 <- log(.expr2)
    .value <- cos(.expr3)
    .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
    .grad[, "x"] <- -(sin(.expr3) * (2 * x/.expr2))
    attr(.value, "gradient") <- .grad
    .value
})

better AD in R?

• basic deriv() can only handle basic functions
• Deriv package
• TMB/RTMB (C++), Chau (2022) (Stan), autodiffr package (Julia)

Rcpp

https://csgillespie.github.io/efficientR/performance.html

pi_rcpp <- cppFunction("
double piSugar(const int N) {
    NumericVector x = runif(N);
    NumericVector y = runif(N);
    NumericVector d = sqrt(x*x + y*y);
    return 4.0 * sum(d < 1.0) / N;
}")

benchmarking against previous example

rbenchmark::benchmark(piR(1e6), pi_rcpp(1e6),
                      columns = c("test", "replications", "elapsed", "relative"))

            test replications elapsed relative
2 pi_rcpp(1e+06)          100   1.876    1.000
1     piR(1e+06)          100   3.280    1.748

more on Rcpp

• lots of examples (gallery)
• not much faster than already-vectorized code
• loops don’t hurt!
• easy to incorporate into packages
• advanced linear algebra via RcppArmadillo, RcppEigen